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The self-healing diffusion Monte Carlo algorithm (SHDMC) [Reboredo, Hood and Kent, Phys. Rev. B 
79, 195117 (2009); Reboredo, ibid. 80, 125110 (2009)] is extended to study the ground and excited states 
of magnetic and periodic systems. The method converges to exact eigenstates as the statistical data collected 
increases if the wave function is sufficiently flexible. It is shown that the wave functions of complex anti- 
symmetric eigen-states can be written as the product of an anti-symmetric real factor and a symmetric phase 
factor. The dimensionality of the nodal surface is dependent on whether phase is a scalar function or not. A 
recursive optimization algorithm is derived from the time evolution of the mixed probability density, which is 
given by an ensemble of electronic configurations (walkers) with complex weight. This complex weight allows 
the amplitude of the fixed-node wave function to move away from the trial wave function phase. This novel 
approach is both a generalization of SHDMC and the fixed-phase approximation [Ortiz, Ceperley and Martin, 
Phys Rev. Lett. 71, 2777 (1993)]. When used recursively it simultaneously improves the node and the phase. 
The algorithm is demonstrated to converge to nearly exact solutions of model systems with periodic boundary 
conditions or applied magnetic fields. The computational cost is proportional to the number of independent 
degrees of freedom of the phase. The method is applied to obtain low-energy excitations of Hamiltonians with 
magnetic field or periodic boundary conditions. The method is used to optimize wave functions with twisted 
boundary conditions, which are included in a many-body Bloch phase. The potential applications of this new 
method to study periodic, magnetic, and complex Hamiltonians are discussed. 
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I. INTRODUCTION 

Following the basic prescriptions of quantum mechanics, 
one could potentially calculate any physical quantity. Finding 
the solutions of the Schrodinger equation, the wave functions 
and the associated energies is all that it is required. However, 
the computational cost of obtaining the solutions of many- 
body problems is well known to increase exponentially with 
the number of particles. Minimizing this exponential cost by 
either improved algorithms or an insightful approximation is 
the central paradigm of condensed matter theory. 

Many physical quantities (observables) are only function- 
als of the probability density: the square of the modulus of 
the wave function. Some other observables, like the ground- 
state energy, are only functionals of the electronic density 1 . 
However, very important quantities, such as the current den- 
sity or excitonic transition matrix elements, depend critically 
on the wave function. 

In confined systems, if the Hamiltonian has time reversal 
symmetry, a real-value wave function is well known to ex- 
ists. However, as soon as periodic boundary conditions are 
introduced in the Hamiltonian or a magnetic field is applied, 
the amplitude of the wave function is, in general, complex. 
If a wave function has a complex amplitude, both its mod- 
ulus $(R) and its complex phase e 1 ^ ( - R ^ can depend on the 
many-body coordinate R = {ri, i"2, • • ■ , rjv c }, where Fj is 
the position of electron j and N e is the number of electrons. 

In addition, for fermions, the many-body wave function 
must change sign when the coordinates of any pair Tj, Tk are 
interchanged in R. In principle, if the wave function is real- 
valued, one only needs to determine the exact surface where 
the wave function is zero and changes sign (the node) to find 



the ground-state energy with diffusion Monte Carlo (DMC) 
methods. Any error in the determination of the node results in 
an overestimation of the ground-state energy^. 

The standard DMC method 5 and improvements 6 related 
to it require, as an input, a trial wave function ^t(R) = 
$T(R)e i ^ (R: ', where both the modulus $ T (R) and </>(R) can 
be chosen to be real. The node is the surface S*t(R) in the 
3N e dimensional many -body space where 4'r(R) = 0. 

The cost of a single DMC step in the standard algorithm is 
polynomial in the number of electrons N e , and can be reduced 
to almost linear if localized orbitals are usecfr— . As a conse- 
quence, the existence of an algorithm that finds the required 
node with polynomial cost in N e has been subject of contro- 
vers y 10,11 . It has been argued that one of the most important 
problems in many-body electronic structure theory is to ac- 
curately find representations of the fermion nodesASi 1 , which 
could help in solving the so-called "fermion sign problem." 

In general, a guess of the node St(R) and the phase </>(R) 
can be obtained from mean field or quantum chemistry meth- 
ods [such as density functional theory (DFT), Hartree-Fock, 
or configuration interaction (CI)]. This initial trial wave func- 
tion is often improved using various methods^— within a 
variational Monte Carlo (VMC) context. This standard ap- 
proach depends on the accidental accuracy of the mean field 
to find the node or the possibility to perform accurate CI calcu- 
lations to pre-select a multideterminant expansion for vpy (R) . 
In addition, it can be claimed that a variational optimization 
of the trial wave function energy or its energy variance only 
improves the nodes indirectly^. 

In the last two years, we have developed a method to cir- 
cumvent the sign problem for the ground and low-energy 
eigenstates of confined system a 18 ! 19 . This method was re- 
cently validated in real molecular systems^. We called the 
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method self-healing diffusion Monte Carlo (SHDMC) since 
the nodes are corrected in a DMC context (as opposed to a 
VMC optimization) and the wave function converges to nearly 
exacW^ or state-of-the-art solutions- , even starting from ran- 
dom. This approach is based on the proof— that by locally 
smoothing the discontinuities in the gradient of the fixed-node 
ground state ^p^(R) at <Sr(R), a new trial wave function 
can be obtained with improved nodes. This proof enables 
an algorithm that systematically moves the nodal surface^ 
SV(R) towards that of an eigenstate. The trial wave func- 
tion is self-corrected within a recursive DMC approach. If the 
form of trial wave function is sufficiently flexible and given 
sufficient statistics, the process leads to an exact eigenstate 
many-body wave function-^—. 

The success of the fixed-node approximation 2 used in the 
standard DMC algorithm for real wave functions is related to 
the quadratic dependence of the error in the fixed-node en- 
ergy with the distance between St(R) and the exact node 4 
5(R). Because the probability density goes to zero quadrat - 
ically at S'(R), errors due to small and short wave-length de- 
partures of <St(R) from S(R) do not propagate far into the 
nodal pocket. Since the DMC energy is dominated by the 
average far from the node, DMC tolerates short wave-length 
departures of S't(R) around 5(R). 

However, if the amplitude of wave function is complex, one 
must also determine its phase </>(R). The ground-state energy 
of complex wave functions can be calculated within the fixed- 
phase approximation 6 of DMC (FPDMC). But any error in the 
phase also results in an overestimation of the ground-state en- 
ergy even if the exact nodes are provided 6 . For complex wave 
functions, moreover, the error in the phase can be more dra- 
matic than the nodal error, since the gradient of phase V0(R) 
is sampled everywhere, and in particular in the regions of large 
probability density (see below and Ref. 0). 

Since (i) periodic or infinite systems are dominant in solid 
state physics, (ii) the ability to calculate complex-valued wave 
functions (with current) is crucial to understanding transport, 
(iii) the response of quantum systems to magnetic fields is key 
for basic understanding of correlated phenomena and even ap- 
plications such quantum computation, (iv) most physical sys- 
tems of interest are not confined, and (v) the error in the phase 
affects the result more than the error in the node, solving "the 
phase problem " is, perhaps, as important as solving the sign 
problem. 

In this paper a method is derived to simultaneously obtain 
not only the node but also the complex amplitude of the trial 
wave function for lower energy eigenstates of Hamiltonians 
with periodic boundary conditions or under applied magnetic 
fields. It is shown that if the phase of the wave function is 
a scalar function, there is a 'special' gauge transformation of 
the many-body Hamiltonian where the wave functions is real. 
These wave functions have nodes that are optimized as in orig- 
inal SHDMC method. If the phase can only be expressed by 
multi-valuate function, the nodal surface may have a reduced 
dimensionality but there is no constraint to update the wave- 
function in SHDMC if the nodes are removed. 

The method is applied and validated in a model system 
studied previously^ 2 * where near-analytical solutions can be 



obtained. The scaling of the cost of this new approach is linear 
in the number of independent degrees of freedom of the phase. 
The method is a generalization of both the "fixed-phase" ap- 
proach 6 and the self-healing DMC algorithms developed to 
circumvent the sign problemi^^ -, The amplitude of the wave 
function is free to adjust to the complex weight of the walkers 
in a recursive approach. 

A study of Refs. [l8l - l20l in reverse chronological order (with 
increasing detail) is recommended before reading this article. 
Studying again the seminal fixed-phase paper by Ortiz, Ceper- 
ley and Martin (OCM) 6 and the importance sampling method 
by Ceperley and Alder- 5 is also highly encouraged. 

The rest of the paper is organized as follows: In Section 
HU the SHDMC and FPDMC methods are generalized and 
blended into a new algorithm that optimizes the complex am- 
plitude (and, if there is one, the node) of the trial wave func- 
tion within a DMC approach. As in the case of SHDMC, the 
trial wave function is adjusted recursively within a generalized 
DMC approach. In Section [HI] the generalization of SHDMC 
is applied to a model Hamiltonian with periodic boundary 
conditions. The results are compared with converged CI re- 
sults for the same model. In Section|IV] the Zeeman splittings 
of the ground and excited states of a model system are cal- 
culated and compared with converged CI results. Section IVl 
describes the results obtained with a realistic Coulomb inter- 
action. Finally, Section[VT]discusses the advantages, perspec- 
tives, and possible applications of these methods for many- 
body problems. 



II. A FREE-AMPLITUDE RECURSIVE DIFFUSION 
MONTE CARLO METHOD 

This section shows how one can obtain an improved trial 
wave function * T (R,r(£ + 1)) = (R|*^ +1 ) by applying 
a smoothing operator D and an evolution operator e ™™ 
(during a small imaginary time r) to the trial wave function 
*r(R, t£) = (R|*t> provided before. The limit r' = It -> 
oo is reached recursively as the iteration index I — > oo. 

Following the seminal ideas of OCM 6 -, ^t(R, t') can be 
written 26 as an explicit product of a complex phase and an 
amplitude * t (R,t') = $ T (R)e i0(R) . OCM chose $ T (R) 
to be symmetric (bosonic like), real, and positive, while the 
phase factor e 1( ^ R ) was antisymmetric for particle exchanges. 
However, the symmetry of a the phase factor is arbitrary: a 
symmetric phase factor can be obtained as 



i0(R) 



* t (R,t') 



^(R,r') 



1/2 



(1) 



since both ^t(R, t') and its complex conjugate change sign 
for particle exchanges. Therefore, any eigenstate can also be 
written as the product of a complex-symmetric phase factor 
e 1 0( R ) (iik e the Jastrow factor) and a real function $y(R) 
where the symmetry of $t (R) depends on whether fermions 
or bosons are considered . 

In this work it is proved (see Subsection III Al l that if the 
phase of a fernionic eigenstate is a scalar function, then 
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$t(R) has the same nodal structure than real functions. Oth- 
erwise <l>y(R) might not be zero except for r-j = Tj. Thus, 
the node of the trial wave function 4 r r(R, r') is given in any 
case by $t(R) but the dimensionality of the nodal surfaces 
depend on the phase. 

The evolution for an additional imaginary time r of 
^^(R, t') is given by 



#r(R, t + t) = e" T ^«* T (R, r) 



$ T (R)e i0(R) 
$T(R,r)e^ (R ). 



(2) 
(3) 
(4) 



Equation (0 includes all the time dependence of the wave 
function in $t(R, t), while the phase (f)(TV) remains fixed 6 .. 

In Eq. (0, e~ T ^ FN is the fixed-node evolution operator, 
which is a function of the fixed-node Hamiltonian operator 
U e FN given by 



H 



FN 



■H+oo lim(9{e- d m [S T (R' Jt) -R]} . (5) 



The second term on the right-hand side of Eq. (0 adds an 
infinite potential 6 at the points R with minimum distance to 
any point on the nodal surface d I „[S , T(R', r') — R] smaller 
than e. The fixed-node Hamiltonian is dependent on I since 
the nodes 5t(R', t') change from one iteration to the next. 

In Eq. 0, the many -body Hamiltonian H is given in atomic 
units by 



(6) 



where Aj = A(r_, ) is a vector potential at point r 3 , V(R) 
includes the electron-electron interaction and any external po- 
tential, and Et is a complex^! energy reference, adjusted to 
normalize the projected wave function, that cancels out any 
phase shift resulting from arbitrary gauge choices for Aj (see 
remarks below). 

Using Eq (0, one can easily obtain 



— [$ T (R,r)] 



H.FN& 



-tUfn 



*t(R, t') 



[£7 L (R,r)-Br]$ T (R,T) 



with 



1^ V?$ t (R,t 



E l (R,t) 



*r(R,r) 
ig|A,+V^(R)| 2 + nR) 



(7) 



(8) 



Vj$r(R, T ) 
$ T (R,r) 



[A J+ V,0(R)] 



Vj • [Aj 



Vj0(R)] 



In Eq. 0, £^(R, t) can be a real constant only if 
$T(R)e 1< ^ R ' > is an eigenstate of Hfn- In general, for an ar- 
bitrary trial wave function, El(R, t) is a complex function of 
R. The real part of El(R, t) is given by the first three terms 
in Eq. (0, while the imaginary contribution is given by the last 
one. OCM's fixed-phase approximation results from consider- 
ing only the real part^i of El (R, r ) . With little effort, one can 
obtain Eq. (3) of OCM's work assuming Iiti[El(TV)} = 0, 
which leads to a continuity-like equation for fluids. 

Note that if 0(R) is held fixed and Im[E L (R)] ^ [see 
Eqs. (0, ©, (0, and (0], then $ r (R, r) must not only 
change its modulus but also must be free to drift away from 
the real values as r increases. 

If at r = an initial distribution of N w walkers /(R, 0) is 
generated to be equal to N w |$t(R) \ 2 , within a generalization 
of the importance sampling algorithm of Ceperley and Alder— 
(see below), /(R, r) should evolve in imaginary time as 



/(R,r) =$ r (R)$ r (R,T). 



(9) 
if 



Clearly /(R, r) can be complex for r > 
Im[E L (R,T)] ^ OfseeEq. ©]. 

Replacing Eq. (0 into Eq. (0, and following a procedure 
almost identical to the one used in Ref. |H one obtains 



a/(R,r) 



Or 



where 



T!»0 



5 £{v|/(R,T)-Vj-[/(R,T)F^} 



[^ i (R)- J B T ]/(R,r), 



VjZn|$ T (R) 



(10) 



(11) 



and £x(R) = El(R,0) is the complex local energy con- 
structed using Eq. (0. To obtain Eq. ( TTOb one must to assume 
that 



Vj$r(R,T) 
*r(R,r) 



$ T (R) : 



(12) 



which implies that unlike the standard DMC algorithm^, there 
is an error in Eq. ( [Tol l when r — > oo if \1/t (R, t) is not an 
eigenstate. This is only an apparent limitation since (i) r at 
first can be made as small as required for Eq. (fT2l to be valid, 
(ii) r can be increased later as the wave function improves and 
converges to an eigenstate, (iii) the limit t' —> oo is reached 
by applying this free-amplitude method recursively (see be- 
low), and (iv) r is already limited to be small in SHDMC with 
correlated sampling so that the weights remain close to 1 (see 
below). 

Although Eq. ( TTOb above for /(R, t) is identical to Eq. 
(1) in Ref. |H it now has a slightly more complex interpreta- 
tion as a stochastic process. Each member of an ensemble 
of systems (walker) undergoes (i) a random diffusion caused 
by the zero-point motion and (ii) drifting by the trial quan- 
tum force In |<E>t(R)| 2 [which depends only on <1>t(R) and 
not on the phase], but in variance with Ref. [H (iii) each 
walker carries a complex phase. In a nonbranching algo- 
rithm, the complex weight of the walkers is multiplied by 
cxp{— [El (TV) — Et] 5t} at every time step. 
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Similar to the case of the "simple" SHDMC algorithm (see 
Refs. [l8U20l for details), the weighted distribution of the walk- 
ers can be written as 



1 Nc 

/(R, r) = lim — £ Wl (k)8 (r - R^ 



(13) 



In Eq. (fT3l . R^ corresponds to the position of the walker i at 
step j of N c equilibrated configurations. The complex weights 
Wf (k) are given by 



,-[E*(k)-E T ]T 



with 



(14) 



(15) 



£=0 



where Et in Eq. ( TT4] > is now a complex energy reference peri- 
odically adjusted so that J2i W/ (k) ~ N c and r is fc(5r (A: is 
a small number of steps and St is a standard DMC time step). 

The trial wave function 4 r j-(R, r' + t) for the next itera- 
tion can be obtained as follows: All wave functions can be 
expanded in a basis as 



*T(R,r') = e J ( R )^A Il (r')<i>„(R 



(16) 



In Eq. ( TToT l. represents a truncated sum, {$„(R)} forms 
a complete orthonormal basis of the antisymmetric Hilbert 
spaced, and is a symmetric Jastrow factor. The A n (r') 

are complex coefficients to be defined [see Eq. (l24ll. Note 
that the expressions^ 



$ T (R) = ± 1 /* T (R, t')^(R, t>), and (17) 
0(R) = ln[* T (R, t')/^* t (R, t')]/{2i) + nn (18) 

allow the computation of all the quantities involved in El (R) 
in terms of gradients and Laplacians of $„(R) and J(R). In 
Eq. (fT8l n is an arbitrary integer that changes the Riemann 
branch of the natural logarithm In, but does not contribute to 
the gradient within a branch. The local energy is thus inde- 
pendent on the choice of n but at the Riemann cuts where, 
sometimes, n has to change to make the phase continuous. 
However, the probability of a walker to touch the Riemann 
cut is, in practice, zero. 

From Eqs. ©, © and (fT3l , one can formally obtain 

§ T (R, t' + t) = e^ (R) /(R, t' + r)/$ T (R) (19) 
= (R|e" T?iFN |^T(T')> • (20) 
The local smoothing operator is defined as 

(R'|£>|R) =<J(R',R) (21) 

= ^e J ( R '>$„(R)C(R)e~' 7(R) - 



Applying Eq. (Bil l to both sides of Eq. (fT9] l, using Eq. ( fT3l , 
and integrating over R, one can easily obtain 

*r(R, t'+t) = (R\De- TfiFN |*t(t')> (22) 

= e J < R )]T(A Il (T' + T))<l>„(R), (23) 

n 

with 

(A n (r'+r)) =lV W/(fc)e- J ^) $;(R ' } 7 (Rl) 

(24) 



where Af = Yli=i e 2J( - R ^ normalizes the Jastrow fac- 
tor. t[R) is the standard time-step correction [Eq. (33) in 
Ref.lm 



7 (R) = 



-1 + v/l + 2|v| 2 r 



with v = 



V$ T (R) 
$ t (R) ' 



(25) 



Note that $r(R) includes a Jastrow factor, thus Eq. 
reduces to the one used in the original "simple" SHDMC al- 
gorithm^ for <f>(R) = 0. 

In addition, as suggested by Umrigar— for the ground-state 
SHDMC algorithm^, correlated sampling can be used also 
for walkers with complex weight. One can sample SX n = 
An(V + t) — A„(t'), which results in 



(Xn(r' + T))=X n (T) + {6\ n ) 



(26) 



{SXr. 



A7^' 

i—l 



[wi(k) - i]7(r|). 



These new A„ (r' + r) [Eq. (l26i ll are used to construct a new 
trial wave function [Eq. (TToT ll recursively within DMC. Equa- 
tion d24l i can be related to the maximum-overlap method used 
for bosonic wave functions^. 

The error of (SX n ) is obtained by sampling 



1 w c 
77^ 



-J(R J ) 



[W/(fc)-l] 7 (R0 



(27) 



The truncation of the expansion of the delta function 
[Eq. ( I2TI 1I is a key ingredient in SHDMC since it decides how 
local is the smoothing operator D and prevents noise to ruin 
the quality of the trial wave function. If the absolute value of 
the error of (SX n ) is larger than |A„(t' + t)|/4, the algorithm 
sets A„ (t' + r) equal to zero [which defines the truncation cri- 
terion used in the sums involved in Eqs. ( [ToT l, (f2TT > and 
<G3]. See Ref.[ll for a detailed theoretical justification of the 
truncation procedure and the algorithm used. Briefly here, the 
coefficients X n are sampled at the end of each sub-block of k 
DMC steps. Statistical data is collected for number of sub- 
blocks M before a wave function update. At first, M is set 
to a small number and increased according to the recipe given 
in Ref. [l9j In short, the algorithm detects automatically the 
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dominance of noise when the projection of two successive sets 
of (SX n ) becomes small and multiplies by a factor larger than 
1 the number of sub-blocks M (see Ref. LL9| for more details). 
As a result, the total number of configurations iV c sampled in- 
creases as the algorithm progresses. Therefore the statistical 
error is reduced, and the number of basis functions retained 
in the expansion increases over time. Thus, the smoothing 
operator D tends to the delta function as M increases, which 
allows the SHDMC method to sample the wave function with 
increasing detail. The first quarter of the data in each block, 
following a wave function update, is discarded. 

In practice, the only difference between this new approach 
and the original SHDMC method is the complex weight and 
the limitation for propagation to small r. Therefore, a non- 
branching algorithm for small r has been used (see Ref. Qjl 
for details). However, there are some formal differences on 
the justification of the convergence of the SHDMC method 
that are discussed in the next subsections. 



A. Gauge transformations and nodal structure of complex 
eigenstates 

The phase </>(R) must be continuous at any point of R 
where ^>t(R) 7^ 0. Otherwise, if 4>(R) is not continuous, 
its gradient in the local energy will introduce an effective infi- 
nite potential at the discontinuity that will force ^t(R) = 0. 
In some cases, however, these discontinuities in the phase are 
not physical 30 and they can be removed by changing the Rie- 
mann sheet index n in Eq. (1181 . As a consequence, the wave 
functions can be split into two classes. In the first class differ- 
ent Riemann sheets of <f>(R) are not connected. In that case, 
one can choose as phase a single sheet of the Riemann surface 
The phase in this class can be a continuous scalar function of 
R for every R. The real wave functions, with constant phase, 
are special case of this class. In the second class, the Rie- 
mann surfaces of <f>(R) for different n are connected at the 
Riemann cuts. Thus a continuous <^>(R) can only be described 
by a multi-valuate function of R. 

Eigenstates with a scalar phase: Since Vj x Vj • A(R) = 
for any scalar function A(R), the magnetic field B = V j x Aj 
is invariant for the gauge transformations^ = Aj + 
Vj[A(R) + c(t)], where A(R) is an arbitrary symmetric 
scalar function of R and c(t) is an arbitrary function of r 
(independent of every r, in R). If ^r(R) is selected to be 
an eigenstate of H for a given gauge and if 0(R) is a scalar 
function , then a change in gauge Aj + Vj<5A(R) could be 
readily compensated in the phase by 

4>(R) = 0(R) - SA(R) + c(r), (28) 

without affecting $t(R) [since Aj and Vj</>(R) always ap- 
pear added in Eq. ([§)]. This property is particularly important, 
since implies that for this class of eigenstates of % there is a 
'special' gauge where the wave function is real. 

Note that if one sets 5A(R) = 0(R) in Eq. (|28]i then 
4>(R) = 0. Therefore, the new phase is a constant that 
can be chosen to be zero. The vector potential in this special 



gauge is a many-body object which includes the gradient of 
the many-body phase of the wave function in a single particle 
gauge. 

The norm $(R) is invariant since the effective potential in 
Eq. (O is invariant using Eq. (l28l . This is expected since the 
expectation value of an arbitrary operator 0(H) must be in- 
dependent of the gauge choice for non-degenerate eigenstates. 
In particular, the nodes, which are given by $(R), are also in- 
variant to gauge transformations. Since the new phase is a 
constant, it can be easily shown that in the special gauge, the 
amplitude of those eigenstates has the same structure as the 
trial wave function used in the fixed-node approximation for 
real wave functions. 

SHDMC self adjusts to an arbitrary gauge change <5A(R) 
because </>(R) is modified recursively by a change in the local 
energy in Eq. $Q of the form 

5E l (R,t)= (29) 
£ I Re [(A, + V^(R)) ■ Vj*A(R)] + -|V^A(R)| 2 | 

-^{^^■V^.ivlMfH)}. 

Eigenstates with a multi-valuate phase: On the other hand 
if 0(R) is not a scalar function then, Vj x Vj • </>(R) 7^ 
and, therefore, V j ■ </>(R) cannot be included in the vector po- 
tential without introducing an artificial many-body magnetic 
field. In that case, as pointed out timely by an anonymous ref- 
eree, there might be nodes only where r 3 = r^. Eigenstates 
with this type of nodes, with reduced dimensionality, can be 
found in states with current, degeneracy or magnetic fields. 

Sumarizing, the norm of the complex wave-functions of the 
eigenstates that have a scalar phase in R has the same struc- 
ture as the real wave function used in the fixed-node approach 
because there is a special gauge transformation where the 
wave function is real valued. This property has a formal im- 
portance since it allows extending theorems developed in the 
context of the fixed-node approximation. Instead, if different 
Riemann sheets of </>(R) are continuously connected, a con- 
tinuous phase cannot be described by a single scalar function. 
Then dimensionality of the nodal surface might be smaller and 
limited to the cases in R where r_, = r^. The nodes are an ob- 
stacle for DMC; SHDMC, however, converges to eigenstates 
regardless of the dimensionality of the nodes (see below). 

B. Remarks on the free-amplitude SHDMC method 

Convergence of SHDMC to eigenstates: In III Al it is shown 
that the wave-function of any fermionic eigenstate can be 
factorized into an anti-symmetric real function $y(R), with 
nodes, and a symmetric phase factor (See Eq. (Q]i). The di- 
mensionality of the nodal surface depends on the phase prop- 
erties. Since the amplitude can not change at the node in 
DMC, nodes are the obstacle to overcame by SHDMC. Con- 
vergence is not affected if the initial trial wave function has no 
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nodes because, in this case, the amplitude and the phase can 
evolve with r everywhere. Indeed, SHDMC can be started 
from a linear combination of real and imaginary parts with dif- 
ferent nodes. As noted by a referee, in this class of functions 
two particles can exchange without crossing a node. However, 
the convergence of SHDMC is not affected in theory and it is 
not affected in practice. But, if the phase is a scalar function, 
nodes will develop as the trial wave function converges to an 
eigenstate. In that limit, a kink at the node should appear~ 
until the exact node is found. 

The convergence of the SHDMC approach when the trial 
wave function approaches an eigenstate and shows nodes, is 
based on the proof 18 that locally smoothing the kinks of the 
fixed-node wave function improves the nodes. This proof 
can be trivially extended to a complex wave function (when 
im[£x(R)] 7^ ), breaking the time evolution into a se- 
quence of pure real evolution followed by an imaginary evo- 
lution. If one assumes i"m[Ex(R)] = 0, the present approach 
reduces to SHDMC method, but using the effective poten- 
tial of the fixed-phase Hamiltonian 6 . Thus the best nodes in 
$t(R, t) for a given phase 0(R) can be obtained by running 
SHDMC in a fixed-phase stage. It is trivial to show that the 
phase, in turn, improves if the imaginary contribution is al- 
lowed to evolve a short time from a trial wave function with 
optimal nodes. In principle, a SHDMC fixed-phase stage can 
be propagated to infinite imaginary time r within a branching 
algorithm. The evolution of the phase, instead, is limited to 
short times [for Eq. ( fT2b to be valid]. In this work, however, 
the real and the imaginary parts of the wave function are al- 
lowed to evolve simultaneously during a short time without 
any observed adverse effect on the accuracy. 

Phase and nodal errors: Note in Eq. ^ that as in the fixed- 
phase approach 6 , an effective potential |Vj</>(R) + Aj\ 2 is 
added to V(R), which depends in turn on V0(R) and Aj. 
Thus small errors in the phase <f>(R) have a global impact [in 
particular far from the node, where |$t(R)| 2 >> 0]. In con- 
trast, small errors in $t(R) only slightly displace the node 
and have smaller impact on the energy- 4 [since the local en- 
ergy is seldom sampled because lim $t(R) < &t(R, t) — > at 
5t(R)]- For eigenstates without nodal surfaces the phase is 
the sole source of error. SHDMC provides a method to correct 
both the phase and the nodal error. 

Complex Et-' Note that when c(r) is changed in Eq. (f28T >, 
it just changes the normalization of ^/ T (R) and, if complex, 
introduces a global phase shift; however, c(r) does not affect 
the local energy or other observables. 

For a random trial wave function 4 f j-(R) and an arbitrary 
choice of gauge for Aj, El(R,t) will have both real and 
imaginary components. The average over the walkers' posi- 
tions will have a real contribution, which affects the norm, and 
a complex contribution, which introduces a global phase shift. 
The average of El (R, t) only contributes to c(r). The corre- 
lated sampling approach is obviously more efficient when the 
distribution of complex weights is centered around 1 since the 
error in the coefficients is minimized [See Eq. [27] and use the 
standard expression for the variance]. Et in Eq. < TT~4T > is thus 
complex. The real part of Et renormalizes the wave func- 
tion (that is, keeps the population of walkers constant). The 



imaginary part of Et removes the global phase shift [the av- 
erage complex contribution of El (R, t) that only contributes 
to c(t)]. For a converged trial wave-function, 7m(-Er) — 0. 

Upper bound properties: The present approach should not 
be considered as a method to estimate the energy of the trial 
wave function but instead as a method to optimize the trial 
wave function before a final FPDMC calculation. However, 
since the real part of El(R) corresponds to the fixed-phase 
approximation, the real part of Et also converges to an up- 
per bound of the ground-state energy. This upper bound can 
be higher than the fixed-phase approximation (if the limit of 
t' — > oo is not reached or the basis {$„(R)} is too small). 
Therefore, a standard FPDMC calculation 6 (with branching) 
must be performed to obtain final values for the ground-state 
energy. 

Known limitations and solutions: The present approach is 
inefficient when the energy of the first excited state E FP of 
the fixed-phase Hamiltonian is too close to the ground-state 
energy E FP , since the coefficient of the excited state com- 
ponent of the trial wave function decays as exp[~T(E FP — 
E FP )]. In that regime, satisfactory results can be obtained as 
follows: begining with a small value for M, keep M fixed for 
a number of iterations until the lower energy excitations decay 
and, then release M, allowing the high energy components of 
the wave function to converge. 

Moreover, the approximate excited-state wave functions 
can be calculated (see Ref. [l9t) and the lowest energy lin- 
ear combination can be determined with correlated function 
Monte Carloi 4 ., VMCi£ or directly in FPDMC using a re- 
stricted basis of low energy states. A final alternative is to 
run this free-amplitude SHDMC method with larger r = kSr 
but using a smaller basis given by a few approximated excited 
states. The excited states can be found as described in the next 
section. 



C. Generalization to excited states 

Earlier estimates of excited state energies in the presence 
of magnetic fields have been made by diagonalizing a matrix 
of correlation functions in imaginary tim o 14 ' 32 . In addition, 
calculations of excited states have been reported with the aux- 
iliary field approach 3 - 3 -. The present algorithm, in contrast, is 
almost identical to the SHDMC metho d 18 ' 19 developed for the 
ground state and lower excitations of real wave functions. The 
only relevant difference with Ref. [HI is the complex weight 
of the walkers. Thus the free-amplitude SHDMC method de- 
scribed above for the ground state can be generalized in a 
straightforward way to study excited states as in Ref. [HI 

Readers are encouraged to follow a detailed theoretical jus- 
tification of the excited state algorithm in Ref. Qjl Here only 
some key steps are described [inparticular, note Eqs. OTb and 
( |32l that were omitted in Ref. [19] and are relevant for a non- 
unitary Jastrow factor]. 

As in the importance sampling algorithm 5 , the generaliza- 
tion given by Eq. ( fTOb requires that 4'r(R, t' + r) be 
zero only at the nodes Sr(R, t') of vffj^R, r'), being free 
to change both its modulus and phase elsewhere. Therefore, 
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t'+t) can develop a projection into any lower energy 
state consistent with 5y(R, t'). To obtain an excited state, 
the wave function 4't(R,t' + t) must be projected in the 
subspace orthogonal to the ground state and any other excited 
state calculated before. In alternative approaches such as cor- 
relation function Monte Carlo 3 ^, the orthogonality of excited 
states is achieved by diagonalizing a generalized eigenvalue 
problem. One could argue that the excited states obtained with 
that approach share nodal error of the ground state. One of the 
advantages of SHDMC is that the diagonalization of a large 
matrix of excitations is avoided, which makes possible the 
consideration of a larger number of degrees of freedom. In ad- 
dition, the nodes of each excitations are found independently. 
But in SHDMC, unless special conditions are satisfied^, one 
must calculate lowerlying energy states before attempting the 
calculation of higher excited states. 

A projector is constructed with approximated expressions 



of the v eigenstates ^(R) 
calculated earlier as 



(R|e J $ M ) = e J ( R )$„(R) 



Pu = 



The operator e J in Eq. (130b is the multiplication by a Jas- 
trow. For a non-unitary e J the set {|$^)} is nonorthogo- 
nal. However, the conjugate (dual) basi s 34,35 that satisfies 
(<t IJi \<b m ) = S^m can be obtained statistically as 



(31) 



with ^ = £„ S-l% since / $;(R)$ u (R)dR = 6 n , m , 

The extension of SHDMC to the next excited l^+i) can be 
thought of as the recursive application of the evolution oper- 
ator e~ kSl ~ n FN > the projector P [Eq d30lll. and a smoothing 
operation D [see Eq. < f21~b l to a trial wave function | v I / T ~j / 1 +1 ) 
[see Eq. d35l)l. This procedure can be derived analytically^ as 
follows: 



lim P e~ rW P|^= +1 ) 

T—>00 



lim P TT ( e -(tT'+kST)up\ |^=o } 

t 

lim PTTfe-^V^* 1 ?) |* £ T t + i) 

I 

lim PTT ( De-^^P) |^=°) (35) 



l*3>+l) 



with 



Replacing e in the infinite product in the second line of 
( 3 °) Eqs. <[35])withe- fc<5T?i in the third line generates the same 
projectori&i£. In turn, we proved 1 - that replacing er &T H with 
a large class of local smoothing operators D has the same ef- 
fect on the nodes. The fixed-node Hamiltonian depends on 
the iteration index I because the trial wave function, the node, 
Et, and the phase are different at every iteration. Finally, the 
norm of the projected function can be fixed by adjusting Et 
in every iteration t. 

For states with inequivalent nodal pockets, special care 
must be taken within the algorithm to avoid systematic errors 
(see Ref.— for additional details about the algorithm). 



it 



lim 



JV C - 



1 ^ Wj(k) q> ra (R^) 

«jy/-^ e -.7(R*)¥£(Rj jT /) 



(32) 



III. CALCULATIONS FOR HAMILTONIANS WITH 
PERIODIC BOUNDARY CONDITIONS 



where ^ T (R,, t') is the trial wave function used to evaluate 
earlier the state /i for t' — > oo. Note that the exponential 
involving — J(R^) moves to the denominator in Eq. d32l as 
compared with Eq. d24l i. Since J(R) is real, the phase of 
(•SpIR) must be conjugated to the phase of (R|$ M ). The co- 
efficients should be sampled during the final FPDMC step 
(i.e., when the final excited energy is sampled). 

The projection of the conjugate function | onto earlier 
conjugate states should also be removed to obtain = 
(^jul-Pj-il where PJ is the transpose of P v . Furthermore, 
statistical errors in can be partially filtered by inverting 
the overlap matrix <n 



(<I» 



= V s^ 1 



(33) 



The scalar products resulting from applying P M in Eq. (f30l > 
are given by 



(34) 



Usually periodic boundary conditions in a supercell with 
dimensions a x , a y and a z are set when studying crystalline 
systems that simulate an infinite solid. By using the Bloch 
Theorem^, the trial wave function at r' = It can be written 
as the product of a many-body phase 3 — times a periodic part2& 
W(R) as 

* T (R, t') = e l &? e kr ^W(R) (36) 



with 



W({ri,' 



^nJ) = U{{r x ,--- ,Tj ■ +a, ■•• ,r N J) 



for any j, where a = a x n x \ + a y n y j + a z n z k, with n M being 
arbitrary integers. U (R), in tunij- can be written as a product 
of a multi-determinant expansion times a Jastrow factor. Each 
orbital entering each determinant in U(R) can be expanded in 
plane waves that satisfy periodic boundary conditions. 

The theory developed in Section [TT] can then be applied to 
periodic systems by setting Aj = 0. Note, however, that since 
U(R) is in general a complex function, the phase entering 
in El(R) [see Eq. [8) must include both the phase of U(R) 
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and the many-body Bloch phase. The resulting wave function 
^^(R-! t' — > oo) corresponds, in general, to a state with cur- 
rent, and its solution can facilitate the calculation of transport 
problems in a many-body context 3 -^. The many-body Bloch 
phase component on the trial wave function is often referred 
to in the literature as twisted boundary condition a 39 ' 40 . 



A. The model periodic system 

Until this new development, the author had promised him- 
self to halt calculations in small model system a 18 ' 19 ' 22 . Those 
small model systems, however, while not very realistic, al- 
low comparisons to be performed with fully converged CI 
calculations. In addition, they can be handled with sym- 
bolic programs like Mathematica, which, while computation- 
ally very slow, are an ideal environment for developing new 
methods and comparing the results with nearly analytical val- 
ues. Therefore, small model systems provide an ideal "work- 
bench" for testing new theories and algorithms. On the other 
hand, no method that fails in the simplest case has hopes of 
succeeding in a realistic calculation involving the more chal- 
lenging Coulomb interaction with a large number of electrons. 
Past experience has shown, in contrast, that earlier SHDMC 
developments tested and developed in small models^ could 
be implemented easily in realistic cases without additional 
complications 2 ^. Indeed, calculations using this method in 
QWALKii reproduced the results obtained for the triplet state 
of He for low magnetic fields^ 2 - starting from a random lin- 
ear combination of determinants constructed with the Hartree- 
Fock solutions without magnetic field. All electron calcula- 
tions of atomic systems with tens of electrons are currently 
under progress and will be published elsewhere^ 3 -. 

DMC calculations with a realistic Coulomb interaction in 
periodic systems require a supercell large enough to prevent 
unphysical image interactions between periodic replicas of the 
electrons from dominating the result. For the purpose of test- 
ing the method, however, a model electron-electron interac- 
tion can be chosen, and the system can be made as small as re- 
quired for numerical convenience. For validating the method, 
the Hamiltonian does not need to be strictly realistic; how- 
ever, one must solve the same Hamiltonian with SHDMC and 
an established benchmark method (CI in this case). 

The model studied in this section is related to the one con- 
sidered in Refs. [l8ll22l and[l9l and consists of two spinless 
electrons in a square of side 1 . However, instead of the hard- 
wall boundary conditions used earlier, periodic boundary con- 
ditions are set. 

Basis expansion: The ground state of the noninteracting 
system is degenerate. Two states with zero total momentum 
can be constructed by placing two electrons with opposite mo- 
menta k = ±7ri or k = ±7rj. The basis chosen to expand 
the wave function is an antisymmetric combination of free- 
particle solutions that satisfy periodic boundary conditions, 
which are plane waves of the form 

e 2ni[(n+l/2)x+my] qj\ 



where |n+l/2| < 6 and \m\ < 5, which results in a two-body 
basis with 1516 functions. 

The confining potential and the interaction potential se- 
lected do not mix the directions 1 and j . They are given by 

V(R) = 47r 2 {cos(27rai)+cos(27ra;2)+cos(27ri/i)+cos(27ry2) 
+ cos[27r(a;2 - a*)] + cos[27r(y 2 - yi )})} . (38) 

The first line of Eq. ( l38l corresponds to an external potential 
applied to electrons 1 and 2. The second line plays the role of 
an interaction potential that depends on the difference between 
the electronic coordinates. 

The Jastrow factor is set to zero to facilitate the analytical 
calculation of the matrix elements of V(R), while the kinetic 
energy is a diagonal matrix. The exact diagonalization of the 
Hamiltonian matrix is the CI result. For % = 1, the energy 
difference 4 ^ between the noninteracting ground state and first 
excited states is 47r 2 . Since the interaction energy in Eq. (l38l ) 
is of the same order of magnitude as the kinetic energy, the 
system is in the correlated regime. 

B. Results and discussion 

Figure Q] shows the logarithm projection Lp(n) = 
/n|(vE'^ / |* T (^T))| of the trial wave function \^S t (£t)) onto 
the n eigenstate of the full CI solution ) as a func- 
tion of the recursive iteration index I. The wave function is 
constrained by the basis to have a many-body Bloch phase 
4>(R) = exp[i(k • i 1 ! + k • r 2 )] with k = 0.97r(i + j) (that is a 
twist angle of 1.87r^d£ both in the i and j directions). 

The initial trial wave function |^r(0)} was chosen inten- 
tionally to be of poor quality to demonstrate the strength 
of the method. The coefficients of |^r(0)) corresponded 
to a linear combination of the first 16 full CI eigenstates: 
|^t(0)) = Yin CnV^n 1 )' wnere the coefficients c„ are com- 
plex numbers of modulus 1 /4 and a random phase. Note that 
the initial trial wave function has no nodes but at the coin- 
cidental points because is a linear combination with random 
phase of different eigenstates of the non interacting Hamilto- 
nian with different nodes. The calculation was run for 200 
walkers with 8t = 0.0004 and r = 0.02. The coefficients A n 
were sampled at the end of each sub-block of k = 50 DMC 
steps. At first the number of sub-blocks M sampled before a 
wave function update was set to 20 and increased according 
to the recipe given in Ref. [l9| and briefly in Section[II] There- 
fore, the statistical error is reduced, and the number of basis 
functions retained in the expansion increases over time. As 
a result both the statistical error and the truncation error di- 
minish, and the wave function continues to improve. The final 
iteration included M = 600 blocks. The total optimization 
run cost w 1.5 x 10 5 DMC steps. 

Figure Q] shows in increasingly lighter shading the results 
Lp(n) corresponding to higher excited states. All the projec- 
tions to the first 16 states start from the same value [— ln(4)] 
by construction. The algorithm, at first, increases the projec- 
tion of the lower energy states at the expense of the higher 
ones (thus Lp(n) approaches zero for low n), while the pro- 
jections with higher n (in lighter gray) become smaller and 
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their Lp(n) is increasingly negative. As the algorithm pro- 
gresses further, the projection on lower-energy excitations 
also starts to decay. Finally, Lp(n) becomes increasingly neg- 
ative for all states except the ground state, which approaches 
zero. 

As the number of recursive iterations I increases, the pro- 
jection onto highly excited states becomes negligible. The val- 
ues obtained for L p (n > 0) are, therefore, dominated by sta- 
tistical noise in the sampling. On the right side of Fig. Q] the 
convergence of the wave function is no longer limited by the 
initial trial wave function but by the statistical noise. Statisti- 
cal noise introduces a projection into higher excited states by 
two mechanisms: (i) the coefficients X n of the trial wave func- 
tion expansion include random noise and (ii) the trial wave 
function develops a projection into excited states because it 
is truncated depending on the relative error of A„, which in 
turn depends on TV ^ 18 ' 19 . Accuracy can be increased only by 
improving the statistics (increasing M and iV c ). 

The residual projection of the trial wave function \^t(^t)} 
for iteration £ on the CI eigenstate l^fi 1 ) is defined as 

Z? p = ln(l-K*f l*r(^)>|). (39) 

The final value for the residual projection for the calcula- 
tion in Fig. Q] is below —7. The value obtained for the 
SHDMC energy is -3 1.842(1 3) as compared with a CI value of 
—31.9486. However, the SHDMC wave function retains only 
70 coefficients in the expansion, whereas the CI has 1516. 
The FPDMC energy obtained with this wave function was 
-32.00(2). 

The results shown in Figure[T]demonstrate that the SHDMC 
method with complex weights is able to correct both the phase 
and the nodal structure of the trial wave function. SHDMC 
converges to the ground-state even starting from a poor quality 
wave function with a random phase. 

C. Many-body band structure 

Common electronic structure methods are based on a 
single-particle picture, and the band structure is given by the 
evolution of the energy as a function of the single-particle 
crystalline momentum. In this case, in contrast, the energy of 
many-body states is a function of the many-body Bloch phase 
e i(£f ' k-rj) or the twist angie^o. 

Figure [2] shows the many-body band structure for the 
ground and first excited states as a function of the global crys- 
talline momentum k = k x \ obtained for the same system stud- 
ied in Fig. Q] The calculations were done using the same pa- 
rameters as in Fig. [T]described above. The trial wave function 
for the ground state with k = started from a linear combina- 
tion of the ground and first excited states of the free-particle 
system with A = Ai = l/\/2- For k ^ 0, the initial trial 
wave function for the ground state was constructed using the 
Bloch part of the converged wave function with smaller |k|. 
The initial trial wave function for the first excited state for 
k = was constructed using a linear combination including 
Ao and Ai orthogonal to the ground state. The trial wave func- 
tions for the first excited states for k ^ were constructed 
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SHDMC iteration index 

FIG. 1: Logarithm of the projection of the trial wave func- 
tion into the lowest 16 eigenstates obtained with CI [Lp(n) = 
Zn|{*J 7 |*T(ft-))|] as a function of the SHDMC iteration index t. 
The results correspond to two electrons in the triplet state with pe- 
riodic boundary conditions and a many-body Bloch phase 4>(TV) = 
exp[i(k ■ ri + k • ra)] with k = 0.97r(T + j) (that is a twist angle of 
1.87r-22ii£). Darker symbols correspond to the projection with lower 
energy CI eigenstates. The initial trial wave function was a linear 
combination of the lowest 16 CI eigenstates with coefficients having 
the same modulus and a random complex phase. 
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FIG. 2: Many-body band structure of a periodic model system with 
two electrons in the triplet state obtained with SHDMC (dots) and 
compared with CI (lines). The figure shows the energy of the ground 
state and the first excited state as a function of the global crystalline 
momentum k (i.e., the many-body Bloch phase or twist angle). 



using the Bloch part of a converged previous calculation with 
the closest value of k and using the projector P to orthogonal- 
ize it with the ground state. CI results are shown with lines for 
validation of the SHDMC results in dots. There is a very good 
agreement between the values obtained with Quantum Monte 
Carlo and CI. In general, however, the Monte Carlo values 
have a higher energy than the CI values. This is due to both 
the error in the complex phase and the nodal error since the 
SHDMC wave function only retains w 70 of the 1516 basis 
functions retained in the CI. The energy difference is reduced 
systematically as the algorithm progresses and more coeffi- 
cients A n are retained in the trial wave function. 
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FIG. 3: SHDMC results (dots) obtained for (a) the first excited state 
and (b) the ground state of a model system compared with CI results 
(lines) as a function of the magnetic field. The system consists of 
two electrons in a two-dimensional square in the triplet state. The 
results shown correspond to the ground and first excited states with 
E symmetry that transform as x + iy. The solution that transforms 
as x — iy can be obtained by changing the sign of B. 



IV. GROUND AND EXCITED STATES WITH APPLIED 
MAGNETIC FIELD 

This section describes the results obtained with the gener- 
alization of SHDMC (described in Section [TTJ> for the ground 
and the first excited state of a model system with an applied 
magnetic field. The results are c ompar ed with CI calculations 
in the same model used in Refs. [l8 



A. Model system with magnetic field 

Briefly, the lower energy eigenstates are found for two 
spinless electrons moving in a two-dimensional square with 
a side length 1 and a repulsive interaction potential of the 
form V(r, r') = 87r 2 7 cos [aTr(x — x')} cos [air(y — y')] with 
a = 1/tt and 7 = 4. The many -body wave function is ex- 
panded in functions $„ (R) that are eigenstates of the nonin- 
teracting system. The basis functions in {$„(R)} are linear 
combinations of functions of the form J\ v sin(m l/ 7rx !V ) with 



m v < 7. Converged CI calculations were performed to ob- 
tain a nearly exact expression of the lower energy states of 
the system 4 f n (R) = 5Zm a m^"i(R)- The matrix elements 
involving the magnetic vector potential A (in the symmetric 
gauge) were calculated analytically using the symbolic pro- 
gram Mathematica and were included in the CI Hamiltonian. 
The Jastrow factor was set to zero in the SHDMC run to facil- 
itate a direct comparison between CI and SHDMC results. 

This paper reports results for the triplet case. In the absence 
of a magnetic field, the triplet ground state is degenerate. Its 
orbital symmetry corresponds to the E symmetry of the D4 
group. One of the solutions with E symmetry transforms as 
x and the other as y. Under an applied magnetic field, the 
time reversal symmetry is broken, and the x and y solutions 
are mixed. Under a magnetic field, the ground state can be 
expanded in a basis of functions that transform as x ± iy. The 
energy of the x — iy solution can be obtained from the energy 
of x + iy by changing B to —B. 



B. Results and discussion 

Figure[3]shows energies of the ground state and first excited 
state of the model system as a function of the magnitude of 
the magnetic field B (the curl of the vector potential A). The 
calculations were run using St = 0.00004 and r = 0.002 
and a total number of DMC steps of 10 5 for each calculated 
point. The calculation for the ground state started from the 
noninteracting ground-state solution as a trial wave function. 
The result obtained for B = for the ground and first excited 
states compared well with the ones obtained with the same 
Hamiltonian in the triplet case reported^ in Table I of Ref . [HI 
Note that in this case, the wave function is complex and the 
coefficients have the freedom to be complex. Thus, in contrast 
with Ref. [HI where a real wave function was enforced, here 
the phase was found within statistical error. 

For B 7^ 0, the time reversal symmetry is broken and so 
is the degeneracy of the x ± iy solutions. For higher (lower) 
magnetic fields, the calculation began by using as the initial 
trial wave function the one obtained previously with a lower 
(higher) magnetic field. 

The excited states were obtained using the method outlined 
in subsection lll Cl and described in detail in Ref. [lj| The lines 
show the CI results for reference. The calculation for the first 
excited state with B = started from a linear combination 
of the ground and first excited state of the noninteracting sys- 
tem orthogonal to the interacting ground state calculated ear- 
lier. The initial trial wave functions of the excited states for 
B were taken from the previous calculations with smaller 
\B\ (keeping the wave function orthogonal to the lower en- 
ergy states with the operator P). Clearly, Fig. [3p shows good 
agreement between SHDMC and CI results for the first ex- 
cited state. 

Table U summarizes the values obtained to construct Fig. [3] 
There is an excellent agreement in the calculations obtained 
for the ground state using SHDMC and CI. The SHDMC en- 
ergy values are, within error bars converged FPDMC results 
indicating that the remaining convergence errors in the basis 
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TABLE I: Comparison of the excitation energies obtained for the 
ground and the first excited state of a model system with two spin- 
less electrons and an applied magnetic field (see Fig. |3)- L rp quanti- 
fies the overlap of the wave functions obtained with CI and SHDMC 
[see Eq. J39H. L var is the variance of the modulus of the walkers' 
weights [see Eq. <40H. 
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343.387 (5) 


343.392(4) 


343.390 


-10.8 


-6.0 


0.8 TT 


344.696 (8) 


344.689(6) 


344.704 


-11.7 


-5.8 


1.6 TT 


347.699 (8) 


347.684(5) 


347.697 


-9.4 


-5.2 



— First Excited State — 



B 


Ei (SHDMC) 


Ei (CI) 


7-1 


Lvar 


-1.6 TT 


394.161 (19) 


394.114 


-7.4 


-5.0 


-0.8 TT 


391.532 (12) 


391.504 


-7.9 


-5.4 


-0.4 TT 


389.744 (12) 


389.741 


-9.1 


-5.7 


-0.2 TT 


388.786 (10) 


388.769 


-9.9 


-5.8 


-0.1 TT 


388.253 (13) 


388.265 


-9.8 


-5.7 


0.0 TT 


388.205 (44) 


387.750 


-5.7 


-5.0 


0.2 TT 


386.697 (17) 


386.694 


-8.9 


-5.7 


0.8 TT 


383.415 (14) 


383.407 


-8.9 


-5.4 


1.6 TT 


379.159 (28) 


379.057 


-7.9 


-4.8 



are small. The agreement is less satisfactory for the excited 
states than in the ground state (using the same computational 
time). It its clear that the residual projections are much larger 
for the excited state than for the ground. 

An independent way to measure the quality of the wave 
function is the logarithm of the variance of the modulus of 
the weights given by 




The variance of the weights does not deteriorate as much as 
the residual projection for excited states, which might signal 
that the differences in the wave functions originate because 
CI and SHDMC minimize different things using a truncated 
basis^ 



V. TEST WITH COULOMB INTERACTIONS 

The calculations with Coulomb interactions were per- 
formed in the same system studied for the ground state in Ref. 



TABLE II: SHDMC and FPDMC energies as a function of an applied 
magnetic field for a model system with two electrons in a triplet state 
in a square box with Coulomb interactions. The quality of the wave 
function is measured by L var [see Eq. J40H, 



State 


B 


SHDMC 


PFDMC 


Lvar 





-1.60 TT 


401.65 (2) 


401.67(4) 


-4.1 





-1.26 TT 


401.80 (3) 




-4.2 





-0.80 TT 


401.92 (3) 


401.87(4) 


-4.2 





-0.40 TT 


403.50 (6) 


402.39(7) 


-3.4 





-0.20 TT 


402.97 (4) 


402.60(5) 


-4.1 





0.00 


402.76 (4) 


402.58(3) 


-4.0 





0.40 TT 


403.23 (2) 


403.20(3) 


-4.6 





0.80 TT 


403.87 (3) 


403.73(3) 


-4.2 





1.26 TT 


404.93 (6) 




-3.7 





1.60 TT 


405.54 (9) 


405.16(4) 


-3.8 


1 


-0.40 TT 


465.37 (10) 




-3.0 


1 


-0.20 TT 


468.55 (7) 




-3.5 


1 


0.00 


454.39 (8) 




-3.4 


1 


0.40 TT 


451.76 (8) 




-3.2 


2 


-0.40 TT 


486.89 (7) 




-3.2 



[l8l and for excited states in Ref. [Hbut now with the additional 
ingredient of an applied magnetic field. The more challenging 
triplet (antisymmetric) state was chosen for this study. 

The calculations were run with the same parameters and 
basis as in Fig. [3] and Table|I]but with a Coulomb interaction 
potential of the form V(r, r') = 207r 2 /|r — r/|. Since the 
average of the Coulomb interaction is much larger than the 
single-particle energy differences, the system is in the highly 
correlated regime. 

Table HIl displays the values obtained for the model system 
with Coulomb interactions for the ground state and some exci- 
tations as a function of the magnetic field. The quality of the 
wave function is characterized by the logarithm of the vari- 
ance of the modulus of weights given by Eq. d40b . Note that 
the variance of the weights increased when Coulomb inter- 
actions are considered when compared with the case of the 
model interaction. This is due to the Coulomb singularity and 
the lack of a Jastrow factor. While the variance of the weights 
is larger in the Coulomb case, the quality of the wave func- 
tion improves from one SHDMC recursive iteration to the next 
(see below). 



A. Improvement of the wave function's node and phase with 
SHDMC 

Figure [4] shows the evolution of the real (a) and the imag- 
inary (b) parts of the local energy El(R) as a function of 
the DMC step for the first excited state of two electrons in a 
square box with an applied magnetic field of OAtt. The cal- 
culations started with a trial wave function with two nonzero 
coefficients chosen to be orthogonal to the ground state cal- 
culated earlier. It can be clearly seen in Fig. |4ja) that as 



i n 




DMC step FIG. 5: (Color online) Logarithm of the variance of the modulus 

of weights as a function of the SHDMC block index I, The results 
correspond to the first excited state with B — 0Att and Coulomb 
interactions (the run shown in Fig.@. 



VI. SUMMARY AND PERSPECTIVES 



FIG. 4: (Color online) (a) Average of the real part of El(R.) ob- 
tained with 200 walkers as a function of the DMC step. The results 
correspond to the first excited state that transforms as x + iy with 
E symmetry of the group D4 of two electrons with Coulomb inter- 
actions in a square box with an applied magnetic field^ B — OAtv. 
(b) Average of the imaginary part of El(R) as a function of the 
DMC step for the same case. The vertical lines mark the end of the 
SHDMC block when the wave function is updated. 



the SHDMC algorithm progresses, the real part of the local 
energy is quickly reduced and stabilized at the first excited 
state energy. The imaginary part of the £x(R) should be 
zero for an eigenstate; otherwise, the divergence of the cur- 
rent is nonzero^. In SHDMC this strong condition is satisfied 
only as the number of recursive iterations, the number of con- 
figuration sampled N c , and the size of the basis Nt, retained 
in the wave function tend to infinity. Figure Etb), however, 
clearly shows that the variance of the raw data obtained for 
im[£x(R)] is reduced as the SHDMC algorithm progresses. 
This is a clear indication of improvement of the phase of the 
wave function. 

Figure [5] shows the evolution of logarithm of the variance 
of the weights [see Eq. d40b l as a function of the SHDMC 
block index (the number of wave function updates). The re- 
duction in weight variance is a clear indication of convergence 
of the trial wave function towards an eigenstate of the Hamil- 
toniani^. 



A method that allows the calculation of the complex am- 
plitude of a many-body wave function has been presented and 
validated with model calculations. An algorithm that finds the 
complex wave function is essential for any study of many- 
body Hamiltonians with periodic boundary conditions or un- 
der external magnetic fields. The method converges to nearly 
analytical results obtained for model systems under applied 
magnetic fields or periodic boundary conditions, with an ac- 
curacy limited only by statistics and the flexibility of the wave 
function sampled. 

It is found that for some eigenstates, the ones where the 
phase is a scalar function of R, there is a special gauge trans- 
formation in which wave function is real. For this class of 
eigenstates the original proof of convergence of SHDMC ap- 
plies. For complex wave functions some fermionic eigen- 
states may not have nodes. In the latter case, as in the case of 
bosons 19 ' 29 , the convergence of SHDMC is not affected since 
the wave function can evolve everywhere. 

This new approach goes beyond both fixed-phase DMC 6 
and SHDMG^— . As in the real wave function version of 
SHDMC, the method is recursive and the propagation to in- 
finite imaginary time is achieved as the number of iterations 
increases. As in FPDMC, the walkers evolve under an effec- 
tive potential that incorporates the gradient of the phase of the 
trial wave function and the vector potential of the magnetic 
field. But in this new algorithm, in contrast to FPDMC, the 
complex amplitude of the wave function is free to adjust both 
its modulus and its phase. After each iteration, the trial wave 
function is improved following a short time evolution of an 
ensemble of walkers. These walkers follow the equation of 
motion of a generalized importance sampling approach. Un- 
like previous attempts, the walkers carry a complex weight re- 
sulting from eliminating the fixed-phase constraint in the time 
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evolution of the mixed probability density. The modulus of 
the weight can be used to calculate real observables, such as 
the energy. The phase of the weight of the walkers is used to 
improve the phase of the trial wave function in the following 
iteration. As in earlier versions of SHDMC, the modulus of 
the weights is also used to improve, simultaneously, the node 
if there is any and the phase of the trial wave function. 

This free amplitude SHDMC method can be used to calcu- 
late not only the ground state but also low energy excitations 19 
within a DMC context. Comparisons with nearly analytical 
results in model systems demonstrate that the new approach 
converges to the many-body wave function of systems with 
applied magnetic fields or with periodic boundary conditions 
for low energy excitations. 

This recursive method finds a solution to "the phase prob- 
lem" and, if there is any, finds the node at the same time. The 
many-body wave function can be used, in principle, to calcu- 
late any observable. However, in very large systems, when 
convergence with the size of the wave function basis cannot 
be fully achieved, a standard fixed-phase calculation should 
be performed as a final step to obtain a more accurate energy. 

Scaling and cost: An analysis of the minimum cost required 
to determine the node and the phase has to take into account 
the number of independent degrees of freedom of the Hilbert 
space. Arguably, no method could scale better than linear in 
the number of independent degrees of freedom of the prob- 
lem studied; otherwise, some degrees of freedom would be 
dependent from each other. A real-space expansion of the 
many-body wave function with fixed resolution Lr is ideal 
for counting independent degrees of freedom. The resolution 
Lr can be connected to the energy cutoff of the excitations in 
a multideterminant expansion. For a complex wave function, 
each point in the many-body space R has two independent 
degrees of freedom (modulus and phase). If the volume of 
a system is proportional to the number of electrons N e , its 

1/3 

size scales as L m aN e (where a is of the order of the 
Bohr radius as). Taking into account the N e ! permutations of 
identical particles, one finds (L/L R ) 3Ne /NJ independent de- 
grees of freedom for each spin channel to determine the phase 
<p(R). Thus, the number of independent degrees of freedom 
scales as exp{N e (3\og(a/ Lr)) + 1)}. The node Sr(R), if 
there is any, requires one less dimension (which, if the nodal 
surface is not too convoluted, could reduce the number of de- 
grees of freedom by only up to a factor (L/ Lr)). Since the 
number of independent degrees of freedom of the phase in- 
creases exponentially with N e , for a fixed resolution Lr one 
cannot find the phase with an algorithm polynomial in N e . 

This generalization of the SHDMC method, though tested 
in small systems, is targeted to be used in large systems. The 
numerical cost of SHDMC scales linearly with the number 
of independent degrees of freedom of the phase per recursive 
step. However, the number of independent degrees of freedom 



(i.e, the size of the basis expansion) should increase exponen- 
tially with the number of electrons N e for a fixed resolution. 
The accuracy of SHDMC is limited by the size of the basis 
sampled, the statistical error, and the number of recursive it- 
erations. The number of recursive steps required increases 
if the product between r and the lowest energy excitations 
is small. The SHDMC method can be used in combination 
with other optimization approaches to accelerate convergence 
in that limit. 

The scaling of the cost of exact diagonalization methods 
such as CI is at least quadratic with the number of degrees of 
freedom. Often a CI calculation is used to preselect a multi- 
determinant expansion to be improved within a VMC context 
before a final FPDMC run. An advantage of SHDMC is that 
it incorporates the Jastrow in the sampling of the coefficients. 
Thus SHDMC might be more efficient than a CI filtering for 
large systems. The linear scaling of SHDMC suggests that it 
could be the method of choice to optimize the wave function 
phase and nodes for calculations in periodic solids. 

The optimization of many-body wave functions with cur- 
rent in periodic boundary conditions is now possible. There- 
fore, the new method can be used as a tool to perform transport 
calculations including many-body effects. The calculation of 
systems with an applied magnetic field is challenging, even 
in the case of small molecules and atoms and particularly so 
when the magnetic field, the many-body interactions, and the 
kinetic energy are of the same order of magnitude^. The cal- 
culations reported in this paper, though in a simple model, 
suggest that the method can be applied to the study of molec- 
ular or atomic systems in that difficult regime. 

Our recent successful application of the ground-state algo- 
rithm for real wave functions 18 to molecular systems^ sup- 
ports the idea that this generalization of SHDMC can also be 
useful for real ab-initio calculations beyond model systems. 
The implementation of the algorithm in state-of-the-art DMC 
codes has been done. Initial results in atomic systems show 
that the many-body wave function improves, which is shown 
by a reduction of the average local energy, the energy variance 
and the variance of the imaginary contribution to the total en- 
ergy. 

Acknowledgments 

The author would like to thank J. McMinis for an introduc- 
tion to the fixed-phase approximation and P. R. C. Kent, and 
G. Ortiz for a critical reading of the manuscript. The author 
also thanks M. Bajdich for sharing all electron calculations 
in atomic systems using this method as supplemental material 
for the referees prior publication. Research sponsored by the 
Materials Sciences & Engineering Division of the Office of 
Basic Energy Sciences U.S. Department of Energy. 



1 P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). 

2 J. B. Anderson, Int. J. Quantum Chem. 15, 109 (1979). 



3 P. J. Reynolds, D. M. Ceperley, B. J. Alder, and W. A. Lester, J. 
Chem. Phys. 77, 5593 (1982). 



14 



B. L. Hammond, W. A. Lester, Jr., and P. J. Reynolds, Monte 
Carlo Methods in Ab Initio Quantum Chemistry (World Scientific, 
Singapore-New Jersey-London-Hong Kong, 1994). 

D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980). 
G. Ortiz, D. M. Ceperley, and R. M. Martin, Phys Rev. Lett. 71, 
2777 (1993). 

A. J. Williamson, R. Q. Hood, and J. C. Grossman, Phys. Rev. 
Lett. 87, 246406 (2001). 

D. Alfe and M. J. Gillan, Phys. Rev. B 70, 161101(R) (2004). 

F. A. Reboredo and A. J. Williamson, Phys. Rev. B 71, 121 105(R) 
(2005). 

D. M. Ceperley, J. Stat. Phys. 63, 1237 (1991). 
M. Troyer and U. J. Wiese, Phys. Rev. Lett. 94, 170201 (2005). 
W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. 
Mod. Phys. 73, 33 (2001). 

G. Ortiz, and D. M. Ceperley Phys. Rev. Lett. 75, 4642 (1995). 
M. D. Jones, G. Ortiz, and D. M. Ceperley, Phys. Rev. E, 55, 6202, 
(1997). 

A. D. Giiclii and C. J. Umrigar, Phys. Rev. B, 72, 045309 (2005); 
A. D. Guclii, G. S.. Jeon, C. J. Umrigar and J. K. Jain, Phys. Rev. 
B 72, 205327 (2005); G. S. Jeon, A. D. Giiclii, C. J. Umrigar, and 
J. K. Jain, Phys. Rev. B 72, 245312, (2005). 

C. J. Umrigar, J. Toulouse, C. Filippi, S. Sorella, and R. G. Hen- 
nig, Phys. Rev. Lett. 98, 110201 (2007). 

A. Liichow, et al., J. Chem. Phys. 126, 144110 (2007). 

F. A. Reboredo, R. Q. Hood, and P. R. C. Kent, Phys. Rev. B 79, 

195117 (2009). 

F. A. Reboredo, Phys. Rev. B 80, 125110 (2009). 

M. Bajdich, M. L. Tiago, R. Q. Hood, P. R. C. Kent, and F. A. 

Reboredo, Phys. Rev. Lett. 104, 193001 (2010). 

^fjv(R-) is not obtained; the new *I't(R.) is sampled directly. 

F. A. Reboredo and P. R. C. Kent, Phys. Rev. B 77, 2451 10 (2008). 

A complex energy reference stabilizes the run for arbitrary gauge 

choices for the vector potential A. 

Equation (2) in OCM work is only strictly valid if Im[EL(R.)] = 
0, otherwise it is the so-called fixed-phase approximation^. 
A complete basis in the symmetric space can be used to calculate 
bosons. Other symmetries of the wave function can be enforced 
with the selection of the basis^. 

Since r' = rl is essentially an iteration index, it is omitted in the 

trial wave function phase and amplitude for clarity. 

C. J. Umrigar, M. P. Nightingale, and K. J. Runge, J. Chem. Phys. 



99, 2865 (1993). 

C. Umrigar (private communication). 
L. Reatto, Phys. Rev. B 26, 130 (1982). 

The position of the cut in the complex plane of the Riemann sur- 
face into sheets is arbitrary. Therefore, the discontinuities in the 
gradient and the effective potential are no physical if they can be 
removed changing n. 

J. D. Jackson, Classical Electrodynamics third edition (John Wil- 
ley & Sons, Inc) New York, (1998). 

D. M. Ceperley and B. Bernu, J. Chem. Phys. 89, 6316 (1988); B. 
Bernu, D. M. Ceperley, and W. A. Lester, Jr., J. Chem. Phys. 93, 
552 (1990). 

W. Purwanto, S. Zhang, and H. Krakauer, J. Chem. Phys. 130, 
094107 (2009). 

K. Hoffman and R. Kunze, Linear Algebra second edition (Pren- 
tice -Hall) New Jersey (1971). 

E. PrugoveCki Quantum Mechanics in Hilbert Space (Academic 
Press) New York (1981). 

N. W. Ashcroft and N. D. Mermin, Solid State Physics (Saunders 
College Publishing Harcourt Brace College Publishers, 1976). 
G. Rajagopal, R. J. Needs, A. James, S. D. Kenny, and W. M. C. 
Foulkes, Phys. Rev B 51, 10591 (1995). 

R. Krcmar, A. Gendiar, M. Mosko, R. Nemefh, P. Vagner, L. Mitas 
PhysicaE 40,1507 (2008). 

C. Lin, F. H. Zong, and D. M. Ceperley, Phys. Rev. E 64, 016702 
(2001). 

The term "twist angle" was introduced^ to avoid confusion with 
other usages of the term "phase". Here the term "many-body 
Bloch phase" is also mentioned since the Bloch theory is well 
understood outside the many-body field. 

L.K. Wagner, M. Bajdich, and L. Mitas, J. Comp. Phys. 228, 3390 
(2009). 

M. D. Jones, G. Ortiz, and D. M. Ceperley, Phys Rev. A 59, 2875 
(1999). 

F. A. Reboredo, M. Bajdich and P. R. C. Kent (work in progrees). 
The energy unit is fi 2 /(2ml 2 ) and the magnetic unit is given by 

e/(c?im (1/2) ). 

There was a small error in the CI calculations reported in Ref.QjJ 
M. D. Jones, G. Ortiz, and D. M. Ceperley, Phys Rev. A 54, 219 
(1996). 



